set more off
clear
cd "C:\Users\amel\Documents\Political Effects of Land-Use Regulation" /*change path as needed*/
log using nhdata.txt, t replace
use nhdata

gen avg00=(gov00+pres00+sen02)/3*100
gen avg16=(gov16+pres16+sen16)/3*100
replace davg0016=davg0016*100
replace lagmedval=lagmedval/100000
replace popratio10=popratio10*100
replace dunem=dunem*100
egen reg0_s=std(reg0)
egen reg2_s=std(reg2)
gen dreg=reg2_s-reg0_s
egen resid0_s=std(resid0)
egen resid2_s=std(resid2)
gen dresid=resid2-resid0
egen ln0697_s=std(ln0697)
egen ln1506_s=std(ln1506)
replace lagdensity=lagdensity/100
replace dens30=dens30/100
replace dmedinc=dmedinc/1000
replace lagmedinc=lagmedinc/1000
gen lagdenssq=lagdensity^2
gen dens30sq=dens30^2
gen dbach_rc=dbach-6.2
gen dgrad_rc=dgrad-3.2

/*appendix figure 14*/
scatter resid2 ln1506_s if zoning15==1, mlabel(town) mlabsize(tiny) mlabpos(0) msize(0) xtitle("Standardized Population Growth, 2006-15") ytitle("Excess Price, 2012-17") msymbol(none) || scatter resid2 ln1506_s if zoning15<1, mlabel(town) mlabsize(vsmall) msize(0) mlabpos(0) mlabcolor(gray) || lfit resid2 ln1506_s if zoning15<1

sort id
sum x_c y_c
spmat idist w2 x_c y_c, id(id) norm(row) vtrunc(0.00000922) replace
spmat export w2 using w2.txt, replace
spmat sum w2, links
/*employment-weighted idm matrix*/
spwmatrix import using idm.csv, wname(W) text connect

splagvar davg , wname(W) wfrom(Stata) ind(reg0_s avg00 resid0_s lagmedval dunem lagdensity ln0697_s lagdenssq dens30 dens30sq island lagmedinc bach_25_2000 grad_25_2000 dbach_rc dgrad_rc bach_25_1970) order(3)

/*is spatial error component required? yes*/
reg davg reg0_s wx_reg0_s avg00 wx_avg00, robust
predict resid, resid
spatgsa resid, w(W) moran geary
drop resid

/*Table 1*/
spreg gs2sls davg reg0_s wx_reg0_s, elmat(w2) id(id) het
spreg gs2sls davg reg0_s wx_reg0_s avg00 wx_avg00, elmat(w2) id(id) het
spreg gs2sls davg reg0_s wx_reg0_s avg00 wx_avg00 lagdensity wx_lagdensity, elmat(w2) id(id) het

/*Table 2, models 5-6*/
spreg gs2sls davg reg0_s wx_reg0_s avg00 wx_avg00 lagdensity wx_lagdensity island, elmat(w2) id(id) het
spreg gs2sls davg reg0_s wx_reg0_s avg00 wx_avg00 lagdensity wx_lagdensity island lagmedinc, elmat(w2) id(id) het

/*excess price instead of supply inelasticity*/
spreg gs2sls davg resid0_s wx_resid0_s avg00 wx_avg00, elmat(w2) id(id) het

/*table 2, model 7*/
spreg gs2sls davg reg0_s wx_reg0_s bach_25_1970 avg00 wx_avg00 dens30 dens30sq wx_dens30 wx_dens30sq, elmat(w2) id(id) het

/*table 2, model 8*/
set seed 244133941

gen davg0= davg - rnormal(-.07292,.0524044)*bach_25_2015
qui spreg gs2sls davg0 reg0_s wx_reg0_s avg00 wx_avg00, elmat(w2) id(id) het
matrix b=e(b)
matrix V=e(V)
matrix S=vecdiag(V)
matmap S E, map(sqrt(@))
matrix rho=e(rho_2sls)
drop davg0

matrix list b
matrix list E
matrix list rho

set matsize 1600

local i=1
while `i' < 200 {
	gen davg`i' = davg - rnormal(-.07292,.0524044)*bach_25_2015
	qui spreg gs2sls davg`i' reg0_s wx_reg0_s avg00 wx_avg00, elmat(w2) id(id) het
	matrix b`i' = e(b)
	matrix b=(b\b`i')
	matrix V`i' = e(V)
	matrix S`i'=vecdiag(V`i')
	matmap S`i' E`i', map(sqrt(@))
	matrix S=(S\S`i')
	matrix E=(E\E`i')
	matrix rho`i' = e(rho_2sls)
	matrix rho=(rho\rho`i')
	drop davg`i'
	local i = `i' + 1
}
svmat b
svmat S
svmat E
rename b1 regsimb
rename b2 wregsimb
rename b3 avgsimb
rename b4 wavgsimb
rename b5 consimb
rename b6 rhosimb
rename S1 regsimv
rename S2 wregsimv
rename S3 avgsimv
rename S4 wavgsimv
rename S5 consimv
rename S6 rhosimv
rename E1 regsims
rename E2 wregsims
rename E3 avgsims
rename E4 wavgsims
rename E5 consims
rename E6 rhosims
sum regsimb
scalar r1=r(mean)
sum wregsimb
scalar r2=r(mean)
sum avgsimb
scalar r3=r(mean)
sum wavgsimb
scalar r4=r(mean)
sum consimb
scalar r5=r(mean)
sum rhosimb
scalar r6=r(mean)
matrix bs=(r1\r2\r3\r4\r5\r6)
sum regsimv-rhosimv
scalar regsimse=sqrt(.2553767+.2553767/200+.2756404)
scalar wregsimse=sqrt(.5595021+.5595021/200+1.473739)
scalar avgsimse=sqrt(.0246199+.0246199/200+.0018112)
scalar wavgsimse=sqrt(.0358745+.0358745/200+.0103509)
scalar consimse=sqrt(1.638235+1.638235/200+25.7993)
scalar rhosimse=sqrt(.0202135+.0202135/200+.0031131)
matrix Simests=(bs,(regsimse\wregsimse\avgsimse\wavgsimse\consimse\rhosimse))
matrix rownames Simests = Reg00 Wreg00 Avg00 Wavg00 _cons _rho
matrix colnames Simests = b se
matrix list Simests

/*college responds to regulation*/
reg dbach_rc reg0_s
hettest
reg dbach_rc reg0_s, robust nocons
predict resid, resid
spatgsa resid, w(W) moran geary
drop resid
reg dbach_rc reg0_s avg00, robust 
reg dbach_rc reg0_s avg00 lagdensity wx_lagdensity, robust 
reg dbach_rc reg0_s avg00 lagdensity wx_lagdensity bach_25_2000, robust 

/*so does median income*/
reg dmedinc reg0_s
hettest
reg dmedinc reg0_s, robust nocons
predict resid, resid
spatgsa resid, w(W) moran geary
drop resid
reg dmedinc reg0_s avg00, robust 
reg dmedinc reg0_s avg00 lagdensity wx_lagdensity, robust 
reg dmedinc reg0_s avg00 lagdensity wx_lagdensity lagmedinc, robust 

log close
